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cn ABSTRACT 

o 

y—i The observations of radio emission from SNR 1987A can be accounted for on the 

^ basis of diffusive shock acceleration of electrons by the supernova blast wave. However, 

Q>^ with this interpretation the observed spectral index implies that the compression ratio 

of the gas subshock is roughly 2.7 rather than the value of 4 expected of a strong 
shock front. We propose that in SNR 1987A, ions also undergo diffusive acceleration 
Q at the shock, a process that is likely to be rapid. Unlike the electron population, the 

^ accelerated ions can have an important effect on the gas dynamics. We calculate this 

^ coupled gas and energetic particle dynamics on the basis of the two-fluid model, in 

which the accelerated ions provide an additional component to the total pressure acting 
on the fluid. By accelerating and possibly heating the upstream plasma, the initially 
strong shock is modified and a weaker subshock with an upstream precursor results. The 
electrons behave as test particles. They are accelerated at the evolving subshock, escape 
downstream, and emit synchrotron radiation in the swept up magnetic field. Two models 
are considered for the surroundings of the progenitor: that of a freely expanding wind 
of number density n oc r~^, and that of a wind confined by a shell of denser material, 
creating a stagnation zone of roughly constant density beyond the standing shock which 
terminates the free wind. We model the observed radio light curves and the relatively 
steep spectrum of SNR 1987A using similar values for the ion acceleration parameters 
to those used in models of cosmic ray acceleration in older SNRs which can also contain 
high Mach number shocks, and find a good fit for the case in which the termination 
shock is located at about 2 X 10^^ m from the progenitor. 

Subject headings: acceleration of particles, shock waves, supernovae : general, supernovae 
: individual : SN1987A, supernova remnants 
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1. Introduction 

The intensity of the radio emission from the rem- 
nant of the supernova SN1987A has increased steadily 
ever since it was detected in July 1990 (Staveley- 
Smith et al. 1992). This behaviour is unlike that 
observed from other radio supernovae - which gen- 
erally exhibit a rapid rise of radio emission followed 
by a steady decline (Weiler et al. 1986) - and has 
stimulated new ideas on the origin of the radiation 
and the interaction of this supernova with its sur- 
roundings (Chevalier 1992a, and Ball &; Kirk 1992a 
- hereafter BK). The latter paper proposed the dif- 
fusive acceleration of electrons and their subsequent 
synchrotron radiation as the mechanism responsible 
for the radio emission, and modelled the light curve 
up to day 1800 by assuming electron injection takes 
place in two distinct components or clumps. This is 
consistent with later observations (Staveley-Smith et 
al. 1993) showing that at about day 2100 the image 
of the remnant consisted of two hot-spots, or clumps 
of emission, superimposed on a diffuse, more or less 
spherical background. However, the model suggests 
that the picture of the supernova shock as a strong, 
unmodified shock front in a gas of adiabatic index 5/3, 
is overly simplistic. The theory of diffusive shock ac- 
celeration makes a definite prediction of the spectral 
index given the compression ratio of a simple unmod- 
ified shock (see Jones and Ellison 1992 and Blandford 
and Eichler 1987 for reviews). On the basis of the 
spectral index of the radio emission, the model fits of 
BK suggest that the electron acceleration takes place 
at a shock of compression ratio ~ 2.7. 

In a recent paper (Kirk, Duffy &: Ball 1994) we ex- 
tended this model by proposing that the shock front 
is not a simple discontinuity travelling out into the 
undisturbed gas of the progenitor's stellar wind, but 
has a structure which is modified by a substantial 
population of accelerated ions. These particles build 
up a precursor to the discontinuity (or 'subshock') in 
the gas flow and lead to a reduction of its compres- 
sion ratio. Electrons have a mean free path which 
is too short to permit diffusion from the subshock 
into the precursor. They are accelerated by the sub- 
shock which, having a lower compression ratio than 
an unmodified shock (and than the overall precursor- 
subshock structure), leads to a relatively steep elec- 
tron spectrum (Bell 1987 and Ellison &: Reynolds 
1991) in agreement with that inferred from the syn- 
chrotron spectrum. If our model is correct, it has 
important implications not only for the future emis- 
sion of SNR 1987A, but also for the theory of other 
radio supernovae (RSNe) as well as for the theory of 
the origin of cosmic rays (in the energy range below 
10^^ eV). The purpose of this paper is to present a 
detailed analysis of the model and compare its pre- 
dictions with the most recent observations of radio 



emission from SNR 1987 A. 

The properties of the medium into which the shock 
front propagates have a strong influence on the ex- 
pected light curve. The model of BK assumes a freely- 
expanding stellar wind around the supernova, which 
implies a number density upstream of the shock front 
proportional to 7""^ (where v is the radius), and (in 
the absence of reconnection or dynamo action) a mag- 
netic field B oc . However, there is substantial 
evidence that the progenitor of SN1987A was, until 
roughly 10'* years ago, surrounded by a dense, slow 
wind typical of a red giant star (for a review see Mc- 
Cray 1993). The transition to the blue-giant phase, 
in which the star found itself upon explosion, was pre- 
sumably accompanied by an increase of the wind ve- 
locity at the stellar surface from the lOkms"^ ap- 
propriate for a red-giant wind, to about 500kms~^ 
(Chevalier &: Fransson 1987). In such a configuration 
the blue-giant wind is expected to inflate a bubble 
inside the relatively dense red-giant wind (Luo &: Mc- 
Cray 1991; Blondin &: Lundquist 1993). Assuming 
spherical symmetry, one would expect the environ- 
ment of the progenitor to have had an inner zone of 
freely-expanding plasma (the blue giant wind) with 
number density n (x. and field 5 oc 7""^ extending 
out to a shock front at ?■ = 7"^, (known as the termi- 
nation shock), surrounded by a stagnation zone with 
roughly constant density and a magnetic field B oc r 
(the shocked blue-giant wind). The outer boundary 
of the stagnation zone is then formed by a contact 
discontinuity separating the shocked blue-giant wind 
from a dense shell of material swept up from the red- 
giant wind. 

The environment around SN1987A shows a marked 
departure from spherical symmetry. Optical observa- 
tions show a ring of material of radius 6.3 x 10^^ m 
(Jakobsen et al. 1991; Panagia et al. 1991) and a 
bipolar nebula structure extending some 2.5 x 10^^ m 
above and below the plane of the ring (e.g.. Chevalier 
1992b). Our hydrodynamical simulations are confined 
to one spatial dimension and so cannot take such de- 
partures from spherical symmetry into account. How- 
ever, even the most recent observations suggest that 
the radio source is still well inside the ring (Staveley- 
Smith et al. 1993) and so this should not yet be a 
severe restriction. Anisotropies in the radio emission 
itself can be attributed to relatively minor enhance- 
ments in the magnetic field strength or in the electron 
injection rate, which leave the overall hydrodynamic 
picture unaffected. 

If the assumption of a spherically symmetric source 
is reasonable for the shock dynamics, at least until the 
supernova blast wave hits the ring, all that remains is 
to fix the location 7"^, of the termination shock front 
between the freely-expanding wind and the shocked 
blue-giant wind. In his model of the radio emission. 
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Chevalier (1992a) attributes the re-emergence of ra- 
dio emission in July 1990 to the arrival of the super- 
nova blast wave at this point, which implies 7"^, ~ 2 
- 3 X 10^^ m, in agreement with estimates based on 
the strength of the blue-giant wind and the size of the 
cavity it has inflated inside the red-giant wind. How- 
ever, this theory predicts a decrease of radio emission 
roughly one year after encounter, in conflict with the 
observed continuing increase. Thus, the precise posi- 
tion of the termination shock, which is is not known 
a priori, cannot be deduced from Chevalier's theory. 
Consequently, we use our model to examine two pos- 
sibilities: (i) that the shock front has not yet reached 
7"^; and (ii) that the shock front passed 7"^, some 1 - 
3 years after explosion. 

The calculations we present fall into two parts: the 
hydrodynamics of the mixture of thermal gas and cos- 
mic rays, and the calculation of the synchrotron emis- 
sion of electrons accelerated at the evolving subshock. 
For the first part we use the two-fluid approach (Dorfi 
1990 and Duffy, Drury &; Volk 1994). Since the su- 
pernova blast wave is currently in its free-expansion 
phase, we assume the fluid velocity behind the gas 
subshock to be constant in time and develop a nu- 
merical procedure for solving the equations governing 
the structure and evolution of the precursor. This 
method is described in detail in §2. Having found the 
hydrodynamic solution, we use it in §3 together with 
the model of time-dependent diffusive shock accelera- 
tion described by BK to construct the electron distri- 
bution as a function of space and time, and to find the 
general form of the synchrotron emission from the ac- 
celerated electrons. In §4 we compare the light curves 
which result from specific models with observed light 
curves. We conclude with a discussion of the impli- 
cations of our results in §5. 

2. The acceleration of ions by a supernova 
shock front 

Cosmic rays (CRs) with energies up to the 'knee' 
of the CR spectrum (at about 10^^ eV), are widely 
thought to be accelerated by the diffusive shock accel- 
eration mechanism operating in supernova remnants. 
Recent work on this theory indicates that the super- 
nova shock front can, whilst expanding into the inter- 
stellar medium, convert between 10% and 30% of the 
total blast wave kinetic energy into CRs, as required 
to sustain their observed intensity (Drury et al. 1989, 
Dorfi 1990, Jones &; Kang 1990, Berezhko et al. 1993). 
With such high conversion efficiencies the reaction 
of the CRs on the gas dynamics cannot be ignored. 
However, most of the energy is converted into CRs 
during the adiabatic or 'Sedov- Taylor' phase of the 
evolution of the remnant, when the shock front has 
already overrun and set into motion a mass compa- 
rable to that ejected initially by the explosion, and is 



steadily decelerating as it expands further. The very 
young remnant of SN1987A is not in this stage of evo- 
lution. The amount of matter overrun by the shock 
front is very small and the ejecta are still expanding 
with essentially undiminished momentum. Therefore, 
because almost all the energy of the explosion is still 
tied up in the kinetic energy of the ejecta, the to- 
tal amount available for acceleration of CRs is small. 
This does not mean the CRs are dynamically unim- 
portant in SNR 1987A. For them to modify the shock 
substantially it suffices that the pressure they exert 
be comparable to the pressure of the post-shock gas. 

We can make a simple estimate of whether or not 
this effect might be important by considering the 
rise of the CR pressure Pc during the test-particle 
phase, i.e., before the CRs become dynamically im- 
portant. Consider a star which explodes into a ra- 
dial, constant-speed stellar wind which has an 7""^ 
density profile. The magnetic field of such a (highly 
conducting) wind has an Archimedian spiral structure 
(Parker 1958) which rapidly becomes azimuthal and 
thereafter drops off as r~^. In the presence of a turbu- 
lent magnetic field, particles can be scattered across 
the shock picking up a small amount of energy each 
time this occurs. The timescale tacc for acceleration 
is proportional to the diffusion coefficient Kp ~ XmtpV 
where A^fp is the mean free path for scattering of a 
particle moving at speed v. Strictly speaking, A^fp 
refers to transport across the azimuthal field. We 
assume that a turbulent magnetic field of the same 
order as the mean field B is present, and estimate 
the effective mean free path for transport perpendic- 
ular to the mean field (Achterberg &: Ball 1994) to 
be the minimum indicated by the quasilinear theory, 
i.e., the particle gyroradius 7"g ocp/B. Therefore, the 
relatively strong magnetic field near the surface of the 
progenitor of SN1987A (Kirk &; Wassmann 1992) im- 
plies very rapid acceleration to high energies (Volk &: 
Biermann 1988). For this rough estimate of Pc we as- 
sume the cosmic rays are ultra-relativistic, neglect the 
effects of adiabatic expansion on their energy and ne- 
glect also the time dependence of the magnetic field at 
the shock front. It is possible to relax these assump- 
tions, but doing so produces no qualitative change 
in the general picture. Assuming that acceleration 
starts at time to, we calculate the upper energy cut- 
off from Pmax = Pmax/^acc B/v. The momentum 
distribution of the CRs is then just that for diffusive 
acceleration at a strong shock, N[p) oc p~^. It follows 
that the temporal dependence of the CR pressure is 
given by 

^(t)4 3 \poJ ^ ^ ' 

where Q{t) is the rate at which particles are injected 
with momentum po and A[t) is the area of the shock 
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front. For a freely-expanding spherical shock, ^ ~ 
and if we assume a constant fraction 77 of the over- 
taken thermal particles is injected, then Q oc rinA is 
constant for an tx oc 7""^ density gradient. The CR 
pressure then varies according to Pc ~ t~^ln(t/to)- 
Whereas the post shock gas pressure, Pq, falls of 
as Pc will initially [t < l.lto) increase very 

rapidly on a timescale much shorter than the expan- 
sion timescale of the shock, Ts/ts, where 7"s(^) is the 
radius of the shock. When t > 1.7to the CR pres- 
sure decreases with time, but it only decreases as 
rapidly as the gas pressure in the asymptotic limit, 
t ^ to- It is therefore plausible that Pc becomes 
the dominant component of the total pressure during 
the early phase of the remnant's evolution, depend- 
ing on the values of the relevant parameters such as 
77 and Kp. Once this happens the test particle esti- 
mates break down and it becomes necessary to solve 
the self-consistent evolution numerically. 

2.1. The Two-Fluid Treatment 

We develop a numerical solution by considering a 
two-fluid system consisting of thermal gas and CRs. 
All the mass is assumed to reside in the gas, but the 
total pressure is the sum of Pq and Pc . In this pic- 
ture there are three processes that can change Pc in 
a frame comoving with the fluid. The first is com- 
pression: the CRs are tied to the fluid through the 
magnetic field and, therefore, change their energy 
when the fluid flow expands or contracts. To cal- 
culate this effect we need the CR equation of state, 
Pc = (7c — 1)-E'c, where Ec is the CR energy den- 
sity and 7c the adiabatic index which decreases from 
5/3 to 4/3 as ultrarelativistic particles begin to dom- 
inate the spectrum. However, without knowledge of 
the spectrum we cannot calculate 7c exactly. A sim- 
ilar problem applies to the second process that can 
modify Pc: the diffusion of CRs off magnetic irreg- 
ularities, because, in general, the diffusion coefficient 
depends on the individual particle energy. Working 
on a macroscopic level with fluid quantities, we re- 
place Kp[p) by an average quantity Kp weighted over 
the (unknown) energetic particle spectrum. The sim- 
plification to a two-fluid model is thus achieved at 
the price of introducing two undetermined closure pa- 
rameters: 7c and Kp. In the present case the ini- 
tial rapid acceleration will quickly lower 7c to 4/3, 
and we use this lower limit for the CR adiabatic in- 
dex. Our calculations therefore underestimate Pc rel- 
ative to Ec, providing a lower limit on the amount 
of shock modification. For Kp we adopt the model 
used in other two-fluid treatments of CR accelera- 
tion (see Drury, Markiewicz and Volk 1989 and Duffy, 
Drury and Volk 1994), namely Kp[r) = /«p(pmax)/4 
and /«p(pmax) = Pmaxc/3g5(7-). This prescription 
results in the required CR conversion efficiencies at 



older SNRs, and the strong X-ray luminosity of young 
SNRs (Dorfi 1990). 

The third process which affects the CR pressure is 
that of injection. Because the two-fluid model does 
not deal with the number density of CR particles, 
but only with their energy density, injection is mod- 
elled by transferring a small fraction e of the avail- 
able mechanical energy at the subshock into CRs. 
The energy flux of injected particles is then given by 
P; = 0.5e [pi{U, - Uif-p2{U, - U2f], where pi,2 
and Ui^2 are the density and velocity, respectively, 
immediately upstream and downstream of the shock. 
This prescription naturally implies a lower injection 
efficiency for weak shocks. The momentum of freshly 
injected CRs, which is needed for the computation 
of the maximum momentum, is taken to be a multi- 
ple, A, of the post-shock gas sound speed Cj times the 
proton mass m, i.e., po = AmCj. In the following, we 
use values for e and A consistent with those needed to 
explain the observed galactic CRs, assuming an SNR 
origin, namely, A = 2 and e = 1 - 5 x 10~^. 
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t/to 

Fig. 1. — The time-dependence of the cosmic ray 
pressure Pc (solid line), the gas pressure Pq (dashed 
line) and the momentum flux P2vl of matter leaving 
the shock (dot-dashed line) with an injection level of 
e = 5xl0~^. AU are normalised to the upstream ram 
pressure of the freely-expanding wind, pqU^ , which 
falls off as 7" ~ ^ . 

In the free-expansion phase of an SNR, there are 
two shock fronts present: an outer or forward shock 
entering the progenitor's wind, and a reverse shock 
moving into the ejecta. The magnetic field imme- 
diately behind the forward shock is just the com- 
pressed field of the wind, which we expect to be 
predominantly toroidal and inversely proportional to 
the shock radius 7"s, at least out to the beginning of 
the stagnation zone. At constant shock speed, there- 
fore, the magnetic field immediately behind the outer 
shock falls off as . On the other hand, the magnetic 
field in the homologously expanding ejecta varies as 
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Fig. 2. — The time dependence of the fluid parameters for a SN shock propagating into a freely-expanding wind, 
for the same parameters used to produce Figure 1. The panels show: (a) the subshock compression ratio, (b) the 
pre- (solid line) and post (dashed line) subshock temperatures, (c) the weighted proton diffusion coefficient Kp in 
units of 7"s?7s, and (d) the maximum CR momentum in units of mc. 
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r/r^(t) _ _ r/r^(t) 

Fig. 3. — The spatial dependence of the fluid parameters at time t = 30to for the case presented in Figures 1 &: 2. 
The panels show: (a) the fluid velocity U, (b) the gas density pc, (c) the gas pressure Pq, and (d) the CR pressure 
in the precursor Pc- 
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B oc 7""^. This is the field encountered by the reverse 
shock, which moves outwards at a speed slightly lower 
than that of the ejecta. Consequently, the field imme- 
diately behind the reverse shock decays as much 
more rapidly than that behind the forward shock. We 
are primarily interested in the synchrotron radiation 
from SNR 1987A, and so confine our attention to 
particles accelerated in the stronger field at at the 
outer shock; all other things being equal, the syn- 
chrotron emission from the reverse shock should not 
be as bright as that from the forward shock. 

To account for the dynamics of the ejecta in the 
free-expansion phase, we assume the flow velocity of 
the gas immediately behind the outer shock is con- 
stant in time. In fact, in spherical geometry, the speed 
of the shocked gas is not exactly equal to the constant 
speed of the ejecta, which constitute a piston driving 
the flow. This is because the pressure behind the 
outer shock drops off roughly as 7""^. If the shocked 
gas were to maintain constant speed, adiabatic expan- 
sion with 7X oc 7""^ would cause a more rapid decline 
(Pg v~^l^ for expansion of a gas of adiabatic index 
5/3). Instead, the shocked gas settles slowly onto the 
piston. We neglect the small velocity difference this 
implies and solve the coupled equations of CR hydro- 
dynamics upstream of the shock. In this particular 
problem, it is important to have good spatial resolu- 
tion of the precursor. We therefore use a formulation 
in which the spatial grid expands along with the shock 
front. Hereafter we use upper case ?7's to denote fluid 
speeds in the rest frame of the explosion's centre, and 
lower case v's for fluid speeds in the rest frame of the 
subshock. Defining ^ = r/rs{t), the coupled two-fluid 
equations in the upstream region can be written in 
the form 
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where U is the bulk plasma flow speed in the rest 
frame of the explosion, and where 
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with Us the speed of the gas subshock. The first two of 
these equations describe the conservation of mass and 



momentum of the gas, with the CR pressure gradient 
included as an additional force; pc is the gas density. 
Equation (4) describes the adiabatic evolution of the 
gas in the precursor region; S is the entropy per unit 
mass. We have assumed there is no dissipation of en- 
ergy into the gas, i.e., that there is no damping of the 
Alfven waves responsible for the scattering of the cos- 
mic rays. Although it is possible to take account of 
this effect, we expect it to be unimportant when, as 
in this case, the direction of the average field is nor- 
mal to the cosmic ray pressure gradient (Jones 1992). 
Equation (5) describes the processes which affect the 
CR fluid. The three terms correspond to compres- 
sion, diffusion and injection, respectively. Finally, in 
order to compute the effective diffusion coefficient Kp, 
we need the magnetic field intensity as a function of 
position and time. Before explosion, the constant ve- 
locity stellar wind contained a toroidal field B oc 1/r, 
up to the beginning of the stagnation zone at 7"^,, after 
which B OL r. On moving into this configuration, the 
shock front forms a precursor, which sets the plasma 
ahead of the shock into motion, compressing both it 
and its frozen-in magnetic field. Assuming this re- 
mains toroidal, we find that B is governed by the 
equation 

^ ' m . (7) 
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Fig. 4. — Evolution of the subshock compression ratio 
for a SN shock propagating into a freely-expanding 
wind, with CR injection levels of e = 0.005 (solid 
line), 0.003 (dashed line) and 0.001 (dot-dashed line). 

The CR pressure gradient at the shock is obtained 
by integrating (5) across the discontinuity at ^ = 1 to 
give 

l+e 
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which is then used as a boundary condition. All that 
remains is to determine the speed of the gas subshock. 
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Fig. 5. — The time dependence of the fluid parameters for a SN shock propagating into a conflned wind. The 
termination shock is at a radius 7"^, = 67U2to- Otherwise, the wind parameters are as used to produce Figure 1. 
The panels show: (a) the subshock compression ratio, (b) the pre- (solid line) and post (dashed line) subshock 
temperatures, (c) the weighted proton diffusion coefficient Kp in units of 7"s?7s, and (d) the maximum CR momentum 
in units of mc. 



This can be obtained from the Rankine-Hugoniot re- 
lations, extended to include the extracted energy flux 
F[. The upstream flow velocity in the rest frame of the 
subshock vi is then given by the quadratic equation 

Vi-Vi {U1-U2) = (9) 

2 pi 

with 7eff = Tg/[1 + e(TG - 1)], where Pi is the gas 
pressure just upstream of the subshock. ^From vi, we 
then obtain the subshock velocity. Us = Ui — vi. The 
details of the numerical scheme used to solve these 
time-dependent equations are given in Appendix A. 

2.2. The Freely- Expanding Wind 

We solve the two-fluid equations with initial con- 
ditions believed to be relevant to the environment 
around the star Sanduleak —69° 202, the progenitor 
of SN1987A (Chevalier &; Fransson 1987, Storey &; 
Manchester 1987). Consider flrst the freely-expanding 
blue supergiant wind moving with a constant speed 
of 500kms"\ Together with a mass-loss rate of 
10~^ Mq yr~^ and a stellar radius of = 3 x 10^° m, 
this determines the proportionality constant of the 
7""^ density proflle. The undisturbed magnetic fleld 
is set to 5 = r^B^/r where 5, = 2 mT (20 G) (see 
Kirk &: Wassmann 1992). When the star exploded it 
ejected about 1 Mq at a speed of some 22, 500kms~^; 



this then is the speed of the downstream medium U2 
that is kept constant in our solutions for the free ex- 
pansion of the remnant. The two parameters control- 
ling injection: e and A (the latter appears only in the 
calculation of the effective diffusion coefficient Kp) are 
set to 5 X 10~^ and 2 respectively. 

The only remaining unknown is the time to at 
which CR injection begins. At very early times, we 
expect the strong magnetic fleld to be very close to 
its unperturbed spiral form. Further out, small fluc- 
tuations in the speed of the wind of the progenitor 
are likely to have produced a more turbulent fleld. In 
addition the energetic particles themselves can gener- 
ate turbulence. We assume that the scattering of CRs 
and injection becomes efficient when the magnitude of 
the fluctuations in the magnetic fleld approaches the 
average fleld magnitude, but the time at which this is 
likely to occur is essentially unknown. We therefore 
take to as a free parameter. However, as we show in 
Appendix A, to introduces the only physical timescale 
into the problem, and the time-dependent evolution 
we compute scales precisely with this time. 

The evolution of Pc, Pg and p2V'^ downstream of 
the subshock (normalised to the upstream ram pres- 
sure of the freely-expanding wind, PqU^ which falls 
off as 7""^), are shown in Figure 1. Pc exhibits a rapid 
initial rise, followed by a decline in absolute terms 
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which is somewhat slower than the ~ fall off of 
Pg and of pqU^ . This behaviour is in accordance 
with our earlier test-particle estimate, and is charac- 
teristic of shock acceleration of CRs in a stellar wind 
cavity. For the parameter values used here the CR 
pressure eventually dominates the thermal pressure 
at t > lOto and the gas subshock starts to weaken. 

This effect is depicted in Figure 2, which shows the 
evolution of the remaining fluid parameters. Panel 
(a) shows that after about lOOto the subshock com- 
pression ratio drops to a level which is comparable to 
that suggested by BK's theory for the radio emission 
of SNR 1987A (i.e., electron acceleration at the gas 
subshock). The remaining panels show the tempera- 
tures of the pre- and post-shock fluids, the weighted 
diffusion coefficient Kp in units of 7"s?7s and the evo- 
lution of the upper CR momentum cut-off Pmax with 
adiabatic losses included. 

Figure 3 shows the spatial dependence of the hy- 
drodynamic solution at time t = 30to- It is notewor- 
thy that the precursor reaches a size which is compa- 
rable to the radius of the subshock. A simple estimate 
based on the equation of motion (3) shows that the 
subshock compression ratio can be significantly mod- 
ified (i.e., Pc can be of the order of pqU^) only if 
the size of the precursor is a significant fraction of 
the radius of the shock. As Pc becomes dynamically 
important the compression ratio between points up- 
stream of the precursor and downstream of the sub- 
shock increases from 4 to values typically greater than 
10 (fig. 3b) since the diffuse, relativistic protons alter 
the overall equation of state. 

The time dependence of the subshock compression 
ratio is sensitive to the injection parameter e. This 
is shown in Figure 4. The higher the level of CR 
injection, the more rapidly the subshock weakens, and 
the lower the value of p at the plateau. 

2.3. The Confined Wind 

The second situation we investigate consists of a 
SN shock expanding into a bubble of freely-expanding 
wind which is terminated at a shock front, outside 
which there is a constant density stagnation zone. In 
the stagnation zone the velocity of the circumstellar 
material decreases as , and the toroidal magnetic 
field increases as 7". In order to simulate this situa- 
tion, we set up an initial density profile which contains 
a sudden transition from n (x. io n = constant 
at a radius 7"^,. The other parameters are as used 
for the freely-expanding wind case discussed above. 
Once again, the time at which injection starts is the 
only physical timescale, provided 7"^, is expressed as a 
multiple of the velocity of the downstream plasma U2 
times to- 

The time dependence of the fluid parameters in 
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Fig. 6. — The effect of the CR injection level on the 
evolution of the subshock compression ratio, for the 
confined wind case. The results for e = 0.005 (solid 
line), 0.003 (dashed line) and 0.001 (dot-dashed line) 
are shown. 



this case is shown in Figure 5. We have chosen to 
place the termination shock at a radius 7"^, = 67U2to- 
The most important difference between the solution 
shown in Figure 5 and the freely-expanding wind solu- 
tion shown in Figure 2 is the decrease of the effective 
diffusion coefficient Kp brought about by the increase 
in B in the stagnation zone. This leads to a character- 
istically different evolution of the subshock compres- 
sion ratio, in which there are two points of inflection 
separated by a more or less pronounced plateau re- 
gion. In the freely-expanding wind case (Figure 2), 
only one point of inflection is seen. 

This can be seen in more detail in Figure 6 which 
shows the effect of the injection parameter e on the 
evolution of the subshock compression ratio. At high 
injection rates, the plateau at a compression ratio of 
about 2.7 is quite persistent. 

The position of the termination shock also has a 
strong effect on the evolution of the subshock com- 
pression ratio, as can be seen in Figure 7. The solid 
line in Figure 7 shows that after a period where the 
subshock compression ratio is essentially constant - 
the 'plateau' discussed earlier - the subshock begins 
to weaken again and eventually disappears. This is a 
general feature of all the solutions for which there is 
substantial CR modification of the shock. The dot- 
dashed line in Figure 7 shows that placing the stag- 
nation zone very close to the star results in the rapid 
smoothing out of the subshock without any significant 
plateau phase. 

3. Acceleration of electrons 

The pressure of the ultrarelativistic electrons re- 
sponsible for the observed synchrotron emission from 
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Fig. 7. — The effect of the the position of the termi- 
nation shock 7"w on the evolution of subshock com- 
pression ratio. The values of 7"^, are 100, 67 and 45 
times U2to (dashed, solid and dot-dashed lines respec- 
tively). 

SNR 1987Ais of the order of 10"® Pa (6x10^ eV cm'^. 
Ball &: Kirk, 1992b), which is only a few percent of 
the cosmic ray pressure. Electrons which are acceler- 
ated by the shock therefore make no contribution to 
the hydrodynamics, and can be treated as test parti- 
cles. Furthermore, the electron diffusion lengthscale, 
at GeV energies, inferred from earlier model fits to the 
observations, Ke/Us ~ 10^^ m (BK), is much smaller 
than the precursor lengthscale which is determined 
by 10'*GeV protons, of ~ 10^^ m obtained in section 
2 (see Figure 3). Therefore electrons in the upstream 
fluid which are overtaken by the shock will initially ex- 
perience the shock precursor as a slow compression, 
then see the gas subshock as a simple, sudden dis- 
continuity in the fluid speed, and subsequently pass 
into the slowly expanding downstream fluid. Diffu- 
sive shock acceleration of electrons will thus occur at 
the gas subshock, which has a much lower compres- 
sion ratio than that of the subshock-precursor system 
responsible for accelerating the CRs. 

Our treatment of the electron acceleration involves 
a straightforward generalisation of that by BK, in- 
cluding now the effect of the evolution of the shock 
strength with time. The two important physical 
processes which modify the electron distribution - ac- 
celeration and adiabatic expansion - are treated sep- 
arately. We assume that the acceleration occurs close 
to the shock on a timescale which is short compared 
to the expansion timescale r/Us- We therefore first 
model the effect of diffusive acceleration in the vicin- 
ity of the shock without adiabatic losses, and then, 
once the accelerated particles escape downstream of 
the shock, the flow is assumed to be diffusion-free and 
adiabatic expansion losses are included. 



In the rest frame of the shock a spatially-averaged 
model of diffusive acceleration leads to the equation 

where, for a particle of speed v, = 4:k[1 / vi + 1 / V2) / v 
is the average time taken to cross and recross the 
shock, A = 4(i)i — V2)/Sv is the average fractional 
increase of particle momentum per cycle, and Pesc = 
4i)2 /v is the probability per cycle of a particle being 
advected downstream away from the shock. Plasma 
flows into the shock at speed vi and out of it at speed 
V2. Equation (10) can be solved by the method of 
characteristics. We assume, for simplicity, that the 
electron diffusion coefficient (and hence tc) is inde- 
pendent of momentum. This is valid for the modelling 
of the radio emission, because the narrowness of the 
"radio window" means that the electrons which con- 
tribute to the radio emission span less than a decade 
in momentum, over which the variation in is rel- 
atively small. We also assume that the electron in- 
jection rate, Q{t), is zero up until some time t^- The 
solution to equation (10) can then be written in inte- 
gral form, the details of which are given in Appendix 
B. When the shock properties are independent of time 
and the injection rate is constant after being switched 
on at ta, the general solution reduces to that of BK, 
namely the general solution reduces to that of BK, 
namely 



N{p,t) 



tcQ 



PoA VPo, 

[h {p -po)-h{p- poe(*-*>)^/*=)] (11) 

where a = Pesc/(2A) = 3/[2(p - 1)] and p is the 
(constant) shock compression ratio. On the other 
hand, when the acceleration takes place at a sub- 
shock which is evolving due to the CR modifications, 
A, tc and Pesc are all functions of time. The gen- 
eral solution to equation (10) can then be integrated 
numerically along the characteristics determined by 
the hydrodynamic solution for the shock evolution. 
The qualitative features of a typical accelerated elec- 
tron distribution resulting from a CR-modified shock, 
which is weakening with time, will then be as fol- 
lows. At momenta just above that of injection, po, 
the spectrum will be a power law as determined by 
acceleration at a simple shock of the appropriate 
instantaneous compression ratio p{t). That is, for 
some range of p just above po, N[p,t) oc where 
q = (p(t) + 2)/(p(t) — 1). At higher momenta the spec- 
trum will be harder than that from a shock of constant 
compression ratio p{t), reflecting the fact that the 
modified shock was stronger at earlier times. Spec- 
tra with just such a concave shape are found when 
equation (10) is integrated numerically using the hy- 
drodynamic solutions of section 2. Figure 8 shows 
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the spectrum which results after 5 years in such a 
case where electron injection is switched on after 2.5 
years and then remains constant, with 1% of electrons 
swept up by the subshock injected into the accelera- 
tion process. The value of the electron diffusion coef- 
ficient is taken to be Ke = lO^^m^s"^. The evolution 
of the hydrodynamics is that of the freely-expanding 
wind solution of Figure 2 with to = 18.3 days. The 
solid line in Figure 8 is the electron spectrum mul- 
tiplied by p^t^y^) with 5(5 yr) = 2.52 corresponding 
to the instantaneous subshock compression ratio of 
p(5yr) = 2.98. The fact that the solid line in Fig- 
ure 8 is essentially horizontal near po indicates that 
the spectrum of the lowest energy electrons is deter- 
mined by the instantaneous subshock compression ra- 
tio. The electron spectrum at the upper momentum 
cut-off is harder, with q = 2.22, and is related to the 
subshock compression ratio at the time when electron 
injection began p(ta) = 3.46 (dashed line). This fol- 
lows from the fact that the electrons with momenta 
near the maximum are those which have spent the 
longest possible time undergoing acceleration, and so 
must have been injected when electron injection be- 
gan. 
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Fig. 8. — Curvature of the electron spectrum pro- 
duced by acceleration at the CR-modified subshock 
of Figure 2, 5 years after the explosion, assuming elec- 
tron injection began at = 2.5 years. The spectrum 
has been multiplied by p^t^^"^) (solid line) to show 
that the most recently injected particles, i.e., those 
near po, have a spectrum determined by the instanta- 
neous subshock compression ratio. The slope of the 
spectrum near the upper cut-off is compared with the 
spectral index given by the subshock compression ra- 
tio at ta (dashed line). 

The second process affecting the electrons is es- 
cape downstream followed by adiabatic losses in the 
expanding fluid. Electrons escape from the shock at 
a rate of N[p, t)Pesc/tc per second, and plasma leaves 
the shock at a speed V2. If the shock is accelerating 



electrons over an area A it follows that the distribu- 
tion function of the escaping electrons is 



47rp^ 



N{p,t) 



AtcV2 
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where the subscript s indicates that this applies only 
immediately downstream of the shock. The diver- 
gence of the downstream flow modifies the distribu- 
tion function which is therefore of the form /(p, r, t) 
and satisfies 



df „^ 1 df 
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where v is the flow velocity of the downstream fluid. 
For radial flow, v = 1/21, this equation may be solved 
using the method of characteristics after introducing 
the Lagrangian (comoving) coordinate R(r, t) which 
is defined implicitly by the equation 

r = R(r,t)+/ dt'v(R,t') (14) 

so that R = r when t = 0. This coordinate is constant 
for a given mass element in the downstream fluid. The 
details of the solution are again given in Appendix B, 
and the general result is 

fip,R,t) 



47r(a!p)2 AtcV2 
N (p[{R+U2t)l{pR)Yl\t) (15) 

where x=[[R+ U2t)/{pR)]^/^ and V2 = Us - U2. 

Since SNR 1987A is still in its free-expansion or 
piston-driven phase, the material behind the shock is 
being pushed out by the ejecta at a speed which is 
essentially constant. There is a slight variation in the 
downstream speed with radius, decreasing marginally 
from the contact discontinuity to the position of the 
shock, but this can reasonably be neglected here, and 
has been neglected in the two-fluid model of section 
2. The relationship between the downstream electron 
distribution /(p, R, t) and the distribution of acceler- 
ated electrons at the shock front N[p,t) is therefore 
unchanged by the inclusion of the cosmic ray mod- 
ifications of the shock front. The only effect of the 
modifications is that on N[p,t) itself. 

Integrating over the downstream electron distribu- 
tion gives their synchrotron flux density at a distance 
D from the source synchrotron flux density at a dis- 
tance D from the source 



/ dpp2 / dR 
Jo JR^it) 

{R+U2tfj,{p)f{p, R,t) 
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(16) 



where Ra,{t) = R[rs{ta,),ta] is the Lagrangian coordi- 
nate of the first injected particles, Rs{t) = R[rs{t),t\ 
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Fig. 9. — The evolution of the pre-subshock magnetic 
field (in mG) at a CR-modified shock propagating 
into a freely-expanding stellar wind (dashed line), and 
into a confined wind which has a termination shock 
at radius 7"^, surrounded by a stagnation zone (solid 
line). 

is the Lagrangian coordinate of the most recently in- 
jected particles, and jv{p) is the single-particle syn- 
chrotron emissivity. There is also a contribution to 
the emitted flux from the accelerating electrons which 
are still in the vicinity of the shock, which is given by 

1 f°° 

To simplify the calculation of the flux density we use 
j^(p) = ao(p/mc)^5^5[i/ — ai[p/mc)'^ B] with ao = 
1.6 X 10-14 W jj^-i rj.-2 andai = 1.3xl0i°Hz T'^. 
The evolution of the flux density with time will de- 
pend on the details of the model but some general 
comments are appropriate here. Emission at a partic- 
ular frequency v will switch on when there are elec- 
trons with Lorentz factor 7l which satisfy the con- 
dition 7L = [i'/aiB{r,t)]'^/^. Therefore, emission will 
appear earlier at low frequencies than at high frequen- 
cies because of the finite acceleration time: emission 
at higher frequencies requires the presence of elec- 
trons with higher energies (or a higher B). After 
switch, on the flux density at a fixed frequency will 
initially increase with time as more electrons are ac- 
celerated up to the required energy. On a somewhat 
longer timescale adiabatic energy losses by the elec- 
trons become important. If the magnetic field (or the 
rate of injection) decreases with radius then emission 
from the newly accelerated electrons will not be able 
to compensate for the adiabatic energy losses of the 
electrons which have already escaped from the shock 
into the downstream plasma. As a result the flux den- 
sity at a given frequency will peak, and then begin to 
fall. 

The evolution of the emission due to acceleration 
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Fig. 10. — Evolution of the rate at which protons 
(and electrons) enter the subshock, for the freely- 
expanding stellar wind case (dashed line), and the 
confined wind case with a termination shock at ra- 
dius 7"w surrounded by a stagnation zone (solid line). 

at a SN shock propagating into, on the one hand, 
an undisturbed freely-expanding stellar wind, and, on 
the other, a confined wind, is quite different. We now 
examine these two situations using the hydrodynamic 
solutions presented in Figure 2 and Figure 5 respec- 
tively. 

The evolution of the magnetic field immediately 
upstream of the subshock, in these two cases, is shown 
in Figure 9. As the shock propagates through the 
undisturbed wind the pre-subshock magnetic field ini- 
tially falls off as 7" - 1 . This is apparent in Figure 9 at 
shock radii 7"s < 0.27"^. When the CR effects become 
dynamically important, the CRs in the precursor be- 
gin to accelerate the stellar wind material ahead of 
the subshock. This compresses the magnetic field, 
and leads to a decline in the pre-subshock field which 
is somewhat slower than 7"-^. In the confined wind, 
the undisturbed field dependence changes at 7"^ from 
B cx to B cx r, but since the magnetic field is 
assumed to be frozen into the plasma, the evolution 
of the field at the subshock is very much smoothed 
out by the modification of the upstream plasma by 
the CR precursor. Thus the confined wind case grad- 
ually evolves away from the freely-expanding wind 
case, starting when the precursor encounters the ter- 
mination shock (at 7"s 0.97"^ in Figure 9) and even- 
tually, when 7"s is considerably larger than 7"^, the 
pre-subshock field increases with increasing 7"s (time). 

The rate at which protons (and electrons) are 
swept up by the subshock is important because we 
assume a certain fraction of these is injected into the 
diffusive acceleration process. The rate goes through 
the same three phases - an unmodified wind phase, a 
wind phase in which CR effects are important, and (if 
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Fig. 11. — Evolution of the flux density from electrons 
accelerated by a CR-modifled SN shock propagating 
into a freely-expanding stellar wind (dashed line), and 
into a conflned wind which has a termination shock 
at radius 7"^, = 2 x 10^^ m surrounded by a stagnation 
zone (solid line). In both cases 1% of electrons inci- 
dent upon the shock are injected into the acceleration 
process. 

applicable) a stagnation zone phase. The evolution of 
the total number of protons (and electrons) crossing 
the subshock per unit time, -ATp^e, is shown in Figure 
10 - again for the two hydrodynamic solutions pre- 
sented in Figures 2 and 5. Initially, the shock expands 
at constant speed, its area increases as 7"^, and since 
the density of the undisturbed wind decreases as 7""^ 
the rate at which plasma crosses the subshock is con- 
stant (7"s < 0.27"w in Figure 10). Then, as CR effects 
become important the plasma upstream of the sub- 
shock is accelerated, and -ATp^e decreases as the shock 
expands further. In the conflned wind case there is 
a gradual turnup in -ATp^e, starting when the precur- 
sor encounters the termination shock [ts ~ 0.97"^, in 
Figure 10). 

The flux density at 843 MHz from electrons ac- 
celerated continuously from a single source region in 
these same two cases - the freely-expanding wind and 
the conflned wind - is shown in Figure 11. We have 
chosen to to be 18.3 days in each case, and for the 
conflned wind we use 7*^ — 67?72^o* Electron injec- 
tion is assumed to begin after 913 days, roughly the 
time that the shock enters the stagnation zone (in 
the conflned wind case), after which 1% of electrons 
crossing the shock are injected into the acceleration 
process. The electron diffusion coefficient is taken to 
be Ke = 10^^ m^s~^. For the case where the shock 
expands into a freely-expanding stellar wind, the flux 
density rises and then slowly decays, for essentially 
the reasons discussed after equation (17). The CR 
effects modify the details of the light curve, but not 
the qualitative trends. However, for acceleration at a 
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Fig. 12. — 'Best' model flt to the observed radio 
emission from SNR 1987 A, for acceleration at a CR- 
modifled shock propagating into a freely-expanding 
stellar wind. For clarity only a subset of the data at 
843 MHz is shown. The uncertainties in the 4.8 GHz 
data are at most the size of the symbols. 

shock expanding into a conflned wind, the flux density 
continues to rise steadily over a very long timescale. 
This reflects the fact that both the magnetic fleld at 
the subshock, and the rate of electron injection (as- 
suming injection of a constant fraction of the elec- 
trons encountered), decrease more slowly after the 
precursor encounters the termination shock than in 
the freely-expanding wind case, and eventually in- 
crease with time. Thus the effect of increased emis- 
sion from freshly-injected electrons continues to domi- 
nate over the decreasing emission from those electrons 
downstream of the shock undergoing adiabatic losses. 

4. Comparison with observations 

To compare the theory with observations we cal- 
culate the electron acceleration and subsequent syn- 
chrotron emission as described in section 3, using the 
two-fluid results of section 2, making now speciflc as- 
sumptions about when and at what rates electrons are 
injected into the acceleration process. Three pieces of 
observational evidence suggest that the continually 
brightening radio emission from SNR 1987A comes 
not from a single source in which the rate of electron 
injection is steadily increasing, but rather that injec- 
tion in discrete components or clumps has switched on 
at different times. Firstly, soon after the remnant was 
detected, the total radio flux density was observed 
to fluctuate on timescales of several days (Staveley- 
Smith et al. 1992), indicating that the emitting re- 
gion was smaller than a spherically-symmetric shell 
of radius equal to that of the shock. Secondly, recent 
high-resolution images of the remnant show that the 
emission itself is clumpy (Staveley-Smith et al. 1993). 
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Fig. 13. — Best model fit to the observed radio emission from SNR 1987A, for acceleration at a CR-modified shock 
propagating into a confined wind. 



Finally there has been at least one marked steepening 
of the light curve (around day 1450) when the rate of 
increase of the flux density jumped rapidly and then 
stayed at its new value, indicating a sudden change of 
injection (Staveley-Smith et al. 1992). We therefore 
model the radio source as two separate components 
in which electron injection starts at different times - 
perhaps on different parts of the shock - and which 
(some time later) produce the switch on or the break 
in the light curves when they commence emitting ra- 
diation at the observed frequency. 

The value of the electron diffusion coefficient is 
constrained by the observed delay in the switch on 
of emission at high frequency (4.8 GHz) compared to 
low frequency (843 MHz). Emission at 4.8 GHz first 
became detectable some 30-60 days after the switch 
on of emission at 843 MHz. ^From equation (7) of 
BK, the ratio of the highest emitted frequencies at 
two successive times ti and ^2 (for acceleration at a 
shock of constant compression ratio) is given by 

^^54 = r^^p[2(ii-i2)A/t.] (18) 

where tc/A is essentially the electron acceleration 
timescale. The observed delay implies that tc/A = 
34 — 68 days. The observed frequency spectra sug- 
gest a value of p = 2.7 for the subshock compression 
ratio, and the two-fluid hydrodynamics solutions of 



section 2 give average values of vi ~ lO^ms"^ for 
the upstream fluid speed around the time that syn- 
chrotron emission was first observed. Together with 
the acceleration time, these values imply a range for 
the electron diffusion coefficient of ~ 1.6 x lO^'^ — 
3 X lO^'^m^s"^. The actual switch on times them- 
selves (as opposed to the delay between them) can 
then be fitted if electron injection begins somewhere 
between 2 years and 2.8 years after the explosion - the 
larger the earlier the start of injection. The val- 
ues of Kg derived here are smaller than those quoted 
in BK because the CR modifications imply a lower 
shock speed for a given ejecta speed. Furthermore, 
the shock continues to weaken and the fluid speeds 
continue to evolve after the switch on of radio emis- 
sion. Therefore, since we assume that the electron dif- 
fusion coefficient doesn't change with time, must 
be somewhat smaller than the above estimates in or- 
der to maintain the required average electron acceler- 
ation time. 

The 'best' fit of the model flux from acceleration 
at a CR-modified supernova shock propagating out 
through a freely-expanding stellar wind is shown in 
Figure 12. The hydrodynamic solution used to pro- 
duce the model light curves in Figure 12 is exactly 
that presented in Figure 2, with to = 18.3 days. The 
electron parameters used are Ke = lO^^m^s"^ for 
both components; tai = 2.5 years and 77ei = 0.008 
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Fig. 14. — Frequency spectra from the 'confined wind' model (lines) versus the observed flux densities (points) at 
the three times for which observed spectra were presented by Staveley-Smith et al. (1992), and at three later times. 
The uncertainties in the data are smaller than the symbols used in the plots. 



for the first component; and ta2 = 2.75 years and 
77e2 = 0.05 for the second. The 843 MHz data plotted 
in Figure 12 (and later figures) are from Ball et al. 
(1994), while the data at 4.8 GHz are from Staveley- 
Smith et al. (1994). 

In this 'freely-expanding wind model' the emission 
from the first source component provides a good fit to 
the data up to the sudden jump in the rate of increase 
of the flux density at around day 1450-1500. Then 
with the addition of a second source component with 
the same value of /«e, an acceptable fit is obtained up 
until about day 1700. This essentially reproduces the 
model fit of BK, and indicates that the CR modifica- 
tions of the expanding shock can explain the relatively 
steep observed spectrum of the radio emission from 
SNR 1987A. However, for the reasons discussed in 
Section 3, the model light curves flatten off and even- 
tually begin to decrease because of the decrease in 
the rate of injection of electrons, and in the magnetic 
field at the subshock. The observed flux however, 
has continued to increase steadily throughout the pe- 



riod of observations to date. We conclude from this 
that electron acceleration at a shock expanding into 
a freely-expanding wind is unlikely to be able to ac- 
count for the continued increase in the radio emission 
from SNR 1987A. Of course, if, for some reason, the 
injection of electrons increases steadily, despite the 
decreasing flux of electrons across the subshock, then 
it would be possible to fit the observations. This case 
is, however, rather artificial - a corresponding contin- 
ual increase in the proton injection rate would lead to 
a departure of the spectrum from that observed. 

Figure 13 shows the 'best' fit of the model flux from 
acceleration at a CR-modified supernova shock prop- 
agating through a confined wind. The hydrodynamic 
solution used to produce the model light curves in 
Figure 13 is exactly that presented in Figure 5 with 
to = 18.3 days and 7"^ = 67U2to- Injection of elec- 
trons from the two source components in this model 
begins after 2.55 years and 2.9 years. The fraction of 
the electrons crossing the subshock that are injected 
are 77ei = 0.006 and 77e2 = 0.02 respectively. The value 
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of Kg is the same in the two source components, and 
is Ke = lO^^m^s"^. This 'confined wind' model pro- 
vides an excellent fit to the observed radio emission, 
clearly reproducing the continuing increase of the flux 
density. The implied radius of the termination shock 
is 2 X lO^^m. In Figure 13 we plot the predicted flux 
up to day 2900, about the time the shock reaches a 
radius of 6.3 x lO^^m, which corresponds to the radius 
of the ring. 

Figure 14 shows the correspondence between the 
frequency spectra obtained from the 'confined wind' 
model, and the observed spectra, at six different 
epochs. The data at 1.5 GHz and 2.4 GHz at the 
three later epochs are from Staveley-Smith (1994). By 
virtue of the reduced compression ratio of the CR- 
modified subshock, the model reproduces and pro- 
vides a physical explanation for the relatively steep 
observed frequency spectra. While the absolute level 
of the radio emission is sensitive to the specific rates 
and timing of the electron injection, the synchrotron 
spectrum is essentially determined by the hydrody- 
namical solution for the CR-modified shock. 

5. Discussion 

The model for the radio emission from SNR 1987A 
due to BK, involving synchrotron emission from elec- 
trons accelerated at the expanding supernova shock, 
provides a good fit to the observations up to about 
day 1800. However it leaves unanswered the ques- 
tion of why the synchrotron spectrum is steeper than 
i/""'^ as expected if the electron acceleration occurs 
at strong shock. In this paper we show that the in- 
clusion of CR effects on the shock structure, via a 
two-fluid treatment, suggests that an initially strong 
SN shock can be substantially modified by the pres- 
sure exerted on the plasma upstream of the shock by 
the CRs which are also accelerated by the shock. This 
process can lead to a relatively rapid evolution of the 
shock. The most important effect for the electron 
acceleration is the formation of a substantial CR pre- 
cursor with a width which is comparable to the shock 
radius. Our estimate of the mean free path for the ra- 
diating electrons is substantially smaller than that of 
the CRs. The electrons therefore see the precursor as 
a gradual change in the fluid properties and undergo 
diffusive acceleration only at the gas subshock. This 
has a much smaller compression ratio than that of 
an unmodified shock, and that of the shock-precursor 
system as a whole, providing a natural explanation 
for the apparent weakness of the shock responsible 
for the electron acceleration in SNR 1987A. 

In order for CR modifications of the SN shock to 
account for the observed spectrum, the efficiency of 
CR injection needs to be sufficiently high for a pre- 
cursor of thickness comparable to the shock radius to 
be established. If the levels of proton injection are too 



low or the value of Kp too high, there is essentially no 
modification of the shock. However, the levels of pro- 
ton injection required to modify the shock sufficiently 
to account for the observed synchrotron spectrum are 
consistent with the levels assumed for CR production 
at supernova remnants in the Sedov- Taylor or adia- 
batic phase of evolution, which are generally believed 
to be the main source of galactic CRs at energies up 
to 10^5 eV. 

The simplest generalisation of the BK model, in- 
cluding CR modification of an initially strong shock 
but retaining the assumption that the shock prop- 
agates through a freely-expanding stellar wind from 
the progenitor star, leads to difficulties in explaining 
the continuing increase in the observed flux densities 
from SNR 1987A. In the original model the flux den- 
sity at a given frequency peaks and then starts to de- 
crease because of the decreasing magnetic field at the 
subshock. When CR modification effects are included 
the decrease occurs sooner, because the CR precur- 
sor accelerates the plasma ahead of the subshock, de- 
creasing the flux of electrons encountering the shock. 
This causes the rate at which electrons are injected 
into the acceleration process to decrease with time, 
assuming injection of a fixed fraction of the electrons 
encountering the shock. The freely-expanding wind 
model (presented in Figure 12) therefore starts to lag 
significantly below the observed flux densities at times 
after about day 1700. 

However, the wind from the blue giant progeni- 
tor of SN1987A may well be confined by an enve- 
lope of more dense material emitted during an ear- 
lier red giant phase. The expanding supernova shock 
may therefore encounter a termination shock in the 
wind, after which the density of the undisturbed cir- 
cumstellar material is independent of radius and the 
magnetic field dependence changes from an de- 
crease to an increase as 7". A model in which the 
electron acceleration occurs at a SN shock expanding 
through such a confined wind, with parameters ap- 
propriate to SNR 1987A, provides a very good fit to 
the observed radio emission at all frequencies. As in 
the original model of BK, two source components are 
needed to account for the sudden jump in the rate 
of increase of the flux density observed around day 
1450-1500. However, an acceptable fit is obtained 
using the same value for the electron diffusion coeffi- 
cient Kg in the the two source components, and with 
models in which the times at which electron injec- 
tion in the two components are within 10% of each 
other, and are close to the time the subshock encoun- 
ters the wind termination shock. The first point is 
important since it suggests that the physical condi- 
tions in the two source components need not be very 
different in order for this model to explain the ob- 
servations. The latter feature implies that the start 
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of electron injection in the two source components 
could occur at the same radius, with the apparent 
time difference in our observer's frame being due to 
light travel time effects. Furthermore, it suggests that 
the start of electron injection may be associated with 
the "impact" of the expanding SN shock on the wind 
termination shock. This encounter is a relatively un- 
spectacular event, even somewhat poorly-defined, be- 
cause the CR precursor of the SN shock smooths out 
the discontinuity in the circumstellar material before 
it is encountered by the subshock itself. Neverthe- 
less, it is tempting to retain the original idea that the 
start of electron acceleration may be the result of the 
change in the conditions at the subshock when it en- 
counters the termination shock of the wind (Ball &: 
Kirk 1992b, Chevalier 1992a). 

The model fits presented here do not unambigu- 
ously determine the position of the wind termination 
shock 7"w, but they do enable us to place limits on 
it. One difficulty which cannot be removed is that of 
the dependence of the results on to, the time when 
CR acceleration begins. However, it is clear that if 
the stagnation zone begins too close to the progen- 
itor, the CR effects will completely smooth out the 
subshock, and hence electron acceleration will cease, 
on a timescale which is so short that it is in conflict 
with the observations of the continuing increase in 
the radio emission. On the other hand, if the wind 
termination shock is too far from the progenitor, the 
electron acceleration at the subshock continues but 
the decreasing magnetic field and electron injection 
again imply that the flux density peaks and starts to 
decline on timescales which are in conflict with the 
continuing increase observed. We therefore conclude 
that the wind termination shock is situated at a ra- 
dius of about 2 X 10^^ m, roughly one third of the 
radius of the dense ring of material around SN1987A. 

Our model for the radio emission suggests that 
the CR pressure behind the shock, i.e., in the down- 
stream plasma, is quite high. When the shock encoun- 
ters the relatively dense ring of material surrounding 
SN1987A, gamma-rays may be produced by the de- 
cay of pions resulting from the collision of CR protons 
with the dense population of protons in the ring. The 
CR pressure is of the same order as the gas ram pres- 
sure (Figure 1) so Pc ~ Piv\. The proton number 
densities of the upstream gas in the precursor and in 
the dense ring are ~ 16cm~^ and ~ 2 x 10'*cm~^ re- 
spectively (Fransson et al. 1989). The rate of gamma- 
ray emission above ITeV can then be estimated from 
Drury, Aharonian and Volk (1994) to be 

N[> ITeV) = q^nEcV ~ 3 x 10^^ s"^ , (19) 

where q~i is a constant, taken to be lO"^'^, contain- 
ing information on the nuclear collision cross sections 
and shape of the CR spectrum and n is the number 



density of protons in the ring which fills a volume 
V. The resulting flux at Earth is 9 x IQ-^^cm-^ s'^ 
somewhat below the threshold of current atmospheric 
Cerenkov detectors. The flux from the confined wind 
can easily be calculated and is slightly lower than this. 
Ideally, a more detailed analysis in which our simula- 
tions are continued until the CR modified precursor 
and subshock hit the dense ring, is required to check 
this estimate. However this is beyond the scope of the 
present paper. It is interesting to note that a similar 
estimate for SN 1993J results in a predicted TeV flux 
which may be detectable (Kirk, Duffy &: Ball 1995). 

The predictive power of our model is limited by 
the known departures from spherical symmetry of 
the distribution of circumstellar material as the SN 
shock approaches the ring and bipolar nebula. Fur- 
thermore, we expect the behaviour of the system to 
change character as the SN shock approaches the ring 
simply because the density of the circumstellar ma- 
terial encountered by the shock will then begin to 
increase rapidly. Our models suggest that the shock 
radius will be comparable to the ring radius some 8 
years after explosion, i.e., in 1995. For the hydrody- 
namic solution of figure 5 in the confined wind case 
the subshock compression ratio will have decreased to 
a value of about 2 at this point. This would imply an 
eventual steepening, though not a disappearance, of 
the radio emission before the point of impact between 
the shock and the circumstellar ring. 
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6. Appendix 
6.1. Appendix A 

All quantities in equations (2) to (7) can be writ- 
ten in terms of their value immediately upstream of 
the shock at to (the time at which CR injection be- 
gins) times a function of r = t/to- For example con- 
sider equation (2) and define p = p/pi[to). With 
fg = r,/r,{to) and A.r] = dr/fg the mass conservation 
equation becomes 



dp 
df] 



[u-iu.) 



dp 



(20) 



This 2-column preprint was prepared with the AAS L^T^X 
macros v3.0. 



where U = U/U[to). Likewise all other equations gov- 
erning the hydrodynamics can be rescaled to include 
only functions of r. 

Equations (1) to (3) describing the dynamics of the 
massive, thermal component of the two fluid system 
are solved explicitly by finite difFerencing. The pre- 
cursor region is divided into M cells each of width 
A^j. If the solution is known at the rath timestep 
then the gas dynamic equations may be difFerenced 
to give the solution a time At later at the [n + l)th 
timestep later at the [n + l)th timestep 



At" 



(.2 jjn + ll2 



Pj+i - Pj ) 



u: 



n + l 
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( pn + l/2 _ pn + 1/2 \ ,„„n 
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At" 



{Sf^^-Sf) (23) 



where xj = ^jU^ — Uj' and the subscript j refers 
to the jth spatial cell. The quantities Uj^^^j^ and 

the solutions of the Riemann problem 
between the jih and (j + l)th cells. Since this scheme 
is explicit it is subject to a Courant Friedrich Lewy 
(CFL) stability limit on the timestep 



At? < 



Hi 



max(C^; + C:-) 



(24) 



with C"j- the local sound speed. Critical to the nu- 
merical accuracy of our method is the resolution of 
the prescursor. Typically we use about 100 cells per 
precursor lengthscale which is sufficient to obtain ac- 
curate results. If we were to use an explicit scheme for 
the CR equation (4) the difFusion term would put a 
strict CFL condition on At that would then scale with 
(A^i)^. Therefore we use a Crank-Nicholson scheme 
where the CR spatial gradient is averaged over the 
timestep. For example the difFusion term in equation 
4 is difFerenced to 



^3- 



2.("+i)^^|A^^ 

-l/2'^j-l/2 ( n + l _ pn + 1 \ 



l/2'^j + l/2 [r}n + l _ pn + A 



+ 
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{Pcj-Pcj-i) 



(25) 



With all other spatial gradients differenced in a sim- 
ilar manner the resulting numerical equation can be 
solved by tridiagonal matrix inversion and provides 
no additional stability constraint on the size of the 
timestep. 

6.2. Appendix B 

The characteristic equation of the acceleration equa- 
tion (10) is dp/dt = pA/tc- Then, for constant and 
Q{t < tg) = 0, the solution to equation (10) can be 
written in the form 



f dt'jiW) 
Jo 

|^*dt'Q(t')exp ^* dt"j{U") 
^exp df'giU") 



(26) 



where g{p,t) = A/t^, h{p,t) = PescAc j{p,t) = 
g{p, t)+p dg{p, t)/dp + h{p, t) and 

^ = P/e^p(_^ dt'g{p,t')^ ■ (27) 

The solution of equation (13) proceeds as follows. 
If we express all functions in terms of the variables 
p, R, t instead of p, r, t then 



(28) 



and defining -D(R, t) = — |V-v we can write equation 
(13) as 



( dfjp, R,t) 
V dt 



R 9P 



The method of characteristics then gives 
ln(p) = / dt' Z)(R,t') + constant 
so introducing 

^ = p/exp (^j^ dt' Z)(R,t') 



p9) 



(30) 



(31) 



it follows that equation (29) can be written in the 
form [df/dt)^ = 0, which has as its solution an ar- 
bitrary function of R and ^(p, R, t). The solution of 
interest follows from the boundary condition which 
relates /[p, R(7" = rs{t),t),t\, where 7"s(i) is the ra- 
dius of the shock, to N[p,t). 



As a specific example we consider the situation dis- 
cussed by BK. where the downstream fluid speed is 
radial and independent of both r and t, i.e. v = ?72 r. 
Equation (14) then reduces to R = r — U2t, and thus 
D{R,t) = -2U2/[S{R+ Uzt)], which on substitution 
into equation (31) gives ^ = p[{R + U2t)/R]^/^. The 
downstream fluid speed U2 is related to the shock 
speed Us by the relation U2 = Us{p— l)/p. Therefore 
if the shock compression ratio is constant it follows 
that the shock speed Us is also constant and thus 
Ts = Ust. The boundary condition is then 



/[p, R{v = Ust,t),t] = fs{p,t) . 



(32) 



The specific solution we require follows from equa- 
tion (12) if we replace p by the combination p[{R + 
U2t)/R]^/^F{R) where F{R) is chosen such that {p[{R+ 
U2t)/R\^'^F{R)}t=R/^u,-u,) = P- Thus F{R) = 
p~^l^ and /(p, R,t) is given by equation (15). 



18 



